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Abstract. The generation and evolution of astrophysical magnetic fields occurs largely 
through the action of turbulence. In many situations, the magnetic field is strong 
enough to influence many important properties of turbulence itself. Numerical simu- 
lation of magnetized turbulence is especially challenging in the astrophysical regime 
because of the high magnetic Reynolds numbers involved, but some aspects of this 
difficulty can be avoided in weakly ionized systems. 

1 Introduction 

The interaction of magnetic fields with turbulence is a basic feature of many 
astrophysical systems. Important, basic MHD turbulence problems common to 
many fields of astrophysics include the nature of MHD turbulence itself, the 
dynamo problem, the effects of magnetic fields on turbulent transport and tur- 
bulent mixing, the effects of small scale turbulence on large scale dynamics, and 
the formation and evolution of current sheets. 

Although progress on all of these problems has been made analytically, nu- 
merical simulations are an increasingly powerful means of approaching them. 
Accurate simulation of astrophysical magnetic fields under turbulent conditions 
presents extreme challenges of its own. The main reason is the smallness of the 
magnetic diffusivity, which leads naturally to the formation of thin current lay- 
ers. Most of the Ohmic dissipation, and most of the change in magnetic topology, 
takes place in these layers. It is important to understand these diffusive effects 
in order to answer questions such as: How do dynamos amplify magnetic fields 
on large scales but not small scales? What determines the ratio of magnetic 
flux to mass in self gravitating regions? How is magnetic energy dissipated in a 
turbulent plasma? 

The purpose of this review is to isolate some important problems related to 
magnetized turbulence in astrophysics, with emphasis on their computational 
aspects. The subject is vast, and we make no claim to be comprehensive. We 
cite literature up to early 2002. 

In §2, we write down the magnetic induction equation, derive some of its 
basic properties, and discuss the use of conservation laws in testing MHD codes. 
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As an illustration, we carry out such tests on the ZEUS code. In §3, we discuss 
the parameter space for astrophysical MHD. In §4, we discuss some results on 
the dynamo problem, the turbulence problem, the formation of singularities, 
and the role of turbulence in large scale dynamics. In §5, we discuss lightly 
ionized media, for which ambipolar drift yields an effective diffusivity which is 
much higher than the Ohmic value. Although the two forms of diffusion are not 
equivalent, ambipolar drift is a partial solution to the problem of small diffusivity. 
Section 6 is a summary and discussion of future prospects. 

2 The Magnetic Induction Equation: Theory and Tests 
2.1 Induction Equation and Consequences 

According to Faraday's law, also known as the magnetic induction equation 

a p> 

- = -cV x E. (1) 

Throughout this paper, we will assume E is given by 

cE = -v x B + a^cJ = -v x B + X n V x B, (2) 

where a is the electrical conductivity, Xq = c 2 /Aira is the magnetic diffusivity, 
and in the last step we have used Ampere's law, neglecting the displacement 
current. 

The first and second terms on the RHS of eqn. (^) represent inductive and 
resistive effects, respectively. In order to compare these terms, we introduce a 
characteristic speed Vq and a characteristic lengthscale Lo, and write v and r in 
terms of dimensionless velocities and coordinates; v = Vqu; r = LqS. Equation 
(||) can then be written as 

cE = -V (u x B - R- 1 V s x B) , (3) 

where the dimensionless parameter R m is defined by 

Rm = -r— ■ (4) 

Xfi 

In some problems it is convenient to set Vq to a typical Alfven speed va = 
B/(4tt / 9) 1 / 2 , in which case R m is known as the Lundquist number and denoted 
by S. 

The case \q oc S^ 1 = is called ideal MHD. It is clear from eqns. (|l|) and 
(^|) that the limit S — > oo is a singular limit, in the sense that at this limit the 
order of the magnetic induction equation drops from second to first. We therefore 
expect that as S — > oo, thin boundary layers, or current sheets, will form. This 



should be anticipated when choosing numerical schemes (§3.3). 

The high S limit describes most astrophysical problems. The ideal form of the 
magnetic induction equation can be solved exactly in terms of fluid trajectories. 
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Suppose that the fluid at position r at time t was at position Tq at time 0. Define 
the deformation matrix D by 

Dij = p-. (5) 
dr 0j 

It can then be shown that the magnetic field B(r,t) is related to the initial field 
Bo(ro(r,t),0) by 

BM)= D ~-"' 0) . (6) 

where \D\ is the determinant of D. Equation (||) is known as the Cauchy solution 
of the magnetic induction equation. It is the basis for the numerical technique 
used by Kinney et al. |53| to simulate 2D MHD turbulence. 

Conservation laws can be used to test MHD codes. The Cauchy solution em- 
bodies magnetic flux conservation. This basic property is in practice difficult to 
test, because in turbulent flow, fluid elements follow complex paths. The Cauchy 
solution only makes sense as long as fluid trajectories do not intersect, but this 
cannot be guaranteed in a finite difference code, which is always somewhat dif- 
fusive. Therefore, we turn to a globally conserved quantity: magnetic helicity. 

The helicity H of a fixed volume V of fluid with magnetic vector potential 
A and magnetic field B = V x A is 

H= [ d 3 rA-B. (7) 
Jv 

Uncurling eqn. (|l|) leads to an evolution equation for A 

F)A 

— =-cE + V<t>, (8) 

where is a free gauge function. According to eqns. (0) and (||), the rate of 
change of H is 

^ = -2c / d 3 rE B+ [ dS ■ (B(f> — cE X A) , (9) 
dt Jv J s 

where we have integrated once by parts and used Gauss's theorem. We now take 
periodic boundary conditions on V , so that the surface integral vanishes, and 
use eqn. (0). Equation (|9|) then becomes 

dH f 

— = -2 / d 3 r\ n B ■ V x B. (10) 
at Jy 

Equation ( |Io[ ) shows that helicity is conserved in an ideal medium. The rate at 
which helicity varies is a global measure of the magnetic diffusivity Af2- Diffusion 
can cause helicity growth as well as decay. 

We obtain a slightly different conservation law if we assume that V is comov- 
ing instead of fixed. In this case, we find that in an ideal medium, H is conserved 
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within V as long as B ■ dS = 0. However, in view of the difficulty of following 
comoving volumes in a turbulent fluid, it is more useful to treat V as fixed, and 
to take it as the computational domain. 

Finally, we derive an equation for magnetic energy W 

w = j yd \^. (11) 

Taking the scalar product of eqn. (|l|) with B, integrating over space, and as- 
suming periodic boundary conditions yields 

dW f J x B / IV x B\ 2 

— = - / d 3 rv / d 3 r\ n l - k i 2 

dt Ur C ./v 4tt 



The first term on the RHS of eqn. ( |12| ) represents the work done by the flow on 
the field, and appears with opposite sign in the evolution equation for kinetic 
energy. The second term represents energy loss by Ohmic decay. 



2.2 Helicity Conservation in the ZEUS Code 



The ZEUS-3D code @,|l00| solves the equations of ideal, compressible MHD 



using a finite difference scheme and a von Neumann artificial viscosity to capture 
shocks. The MHD induction equation is followed using the method of consistent 
transport along characteristics j46). ZEUS-3D is publically available, and has 
been of great service in the astrophysical community. 

If there were no numerical dissipation in the ZEUS code, magnetic helicity 



would be strictly conserved (see eqn. (10)). Here, we investigate helicity conser 



vation in the ZEUS code in two applications. 

The first problem is the evolution of a twisted magnetic flux tube which is 
unstable to the kink mode. Initially, the helicity of the tube is entirely in the 
twist of the field about the axis. Theory predicts that the instability drives the 
system to a new equilibrium state in which some of the helicity is carried by a 
writhing deformation of the tube axis, and the total magnetic energy is reduced 
while the total helicity is fixed. 

The simulations confirm the broad outline of this picture. The growth rate 
of the kink during the linear phase agrees with an analytical calculation, and 
significant motion occurs only during the kinking phase, while magnetic energy 
is being released. However, both magnetic energy and magnetic helicity decline 
steadily once the system reaches equilibrium, as shown in Figure (Q) for com- 
putations with 80 3 and 160 3 gridpoints. During the dynamical kink phase, the 
energy declines much faster than the helicity, confirming the importance of both 
dynamical and resistive processes in the evolution of magnetic energy. 



We have used eqn. (10) to estimate the mean pseudo- magnetic diffusivity 
(Aq), i.e. the mean numerical magnetic diffusivity, by writing the integral on the 
right hand side as 

d 3 r\ n B V x B = (An) / d 3 rB ■ V x B, (13) 
v Jv 
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Fig. 1. Kinetic and magnetic energy, and helicity against time for the twisted magnetic 
flux tube, t is given in units of the Alfven crossing time for the flux tube with diameter 
a. Helicity is given in units of Lq<P 2 /2n, where Lq/2iv is the number of rotations by 2ir 
of the twisted field lines about the tube axis over the domain length L, and $ is the 
total axial flux of the tube. 



which can be regarded as the definition of (A^). The result is shown in Figure 
(U), where (Xn) is plotted for the two simulations in units of Axvao, where Ax is 
the grid scale. The near coincidence of the two curves shows that the numerical 
diffusion is linear in Ax. During the dynamical phase, (Xn) is enhanced over its 
"quiet" value cqAxvao by about a factor of 3, perhaps suggesting that numerical 
resistivity, like numerical viscosity, is proportional to velocity. 

The results shown in Figure (^|) allow us to estimate the Lundquist number S. 
According to eqn. S can be written in terms of eoj the number of gridpoints 
N, and the magnetic lengthscale Lb in units of the box size L as 

S=±f — . (14) 

In these computations, the tube has an initially Gaussian magnetic profile with 
Lb/L = 0.1. Taking eo ~ 0.01, which is representative of the "quiet" value 
shown in Figure (0) and N 1 / 3 - 100, we see that S < 10 3 . 

We have also computed the evolution of H in models of molecular clouds 
[^7|^], with initially uniform magnetic fields, stirred by supersonic turbulence. 
In these models H is initially zero, and should remain so, although the helicity 
density A ■ B can vary arbitrarily between positive and negative values from 
point to point. Figure (||) shows H relative to E m L for two different runs at 
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Fig. 2. Mean magnetic diffusivity against time for the twisted magnetic flux tube. 



various times, measured in units of the sound crossing time. In these units, the 
errors in H are at the 1% level, and, reassuringly, there is no evidence that H 
drifts steadily away from zero. The variance of H around zero is smaller for the 
run at Alfven Mach number Ma = 0.7 than for the run at Ma = 6.7, presumably 
because the field becomes less tangled if it is relatively strong. 

Figure (^) shows mass density and helicity density, for the weak field model, 
in a plane perpendicular to the mean magnetic field. The figure shows that 
helicity density varies strongly with position, and that it is well correlated with 
mass density (the correlation is not as good in the strong field case). 

The density compressions in these models are associated with shock fronts. 
Helicity density can also increase across a shock. For example, consider a locally 
sheared, force free magnetic field of the form 

B = Bq (xcoskoz + ysmkoz) . (15) 

The vector potential is 

A = -k^B, (16) 

so the helicity density is —k^ ) 1 B^ j . Suppose the fluid is compressed in the z 
direction, so that the initial and final coordinates zq and z are related by z — azo, 
with a < 1. It can be shown from eqn. (|^), or from intuitive arguments, that 
B still has the form (|l5|), but with Bq — > o.~~ 1 Bq and fco — » a _1 fco- The helicity 
density then increases by a factor a -1 , which is the same factor by which the 
mass density increases. 
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Fig. 3. Total magnetic helicity against time for two models of turbulent molecular 
clouds at M A = 0.7 and M A = 6.7. 



density helicity 




Fig. 4. Density (left panel) and helicity (right panel) for the weak field model with 
M A = 6.7. 

3 Parameter Space 

The dimensionless parameters of astro-physical turbulence define the scope of the 
problem and should influence the choice of numerical technique. Here we define 
and give quantitative expressions for some important parameters. 

3.1 Macroscopic Parameters 

The ratio of gas pressure to magnetic pressure is usually denoted by (3 = 8ttP/B 2 , 
which is closely related to the ratio of the sound speed vs to the Alfven speed 
va; Vg/v\ = 7/3/2, where 7 is the ratio of specific heats. If (i is either very 
small or very large, the acoustic and Alfven frequencies are very different, with 
consequences for choice of simulation technique. 
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The ratios of the mean flow speed v to vs and va , respectively, are denoted by 
the sonic and Alfven Mach numbers M and Ma- When M ^> 1, compressibility 
effects are important and the flow is pervaded by shocks, requiring an algorithm 
which treats them accurately. 

The importance of self gravity is measured by the Jeans length Aj. Turbulent 
pressure effectively increases A,/. This is captured by a heuristic expression for 
A,/ in a turbulent medium 



TTV 2 S (1+M 2 (\j)) 



Gp 



1/2 



(17) 



where M 2 (Aj) is the Mach number of the turbulence at wavelengths less that 
Aj; turbulence at longer wavelengths does not contribute to pressure support 
[ p2p^ , |9C| . Equation ( |l7| ) reverts to the usual definition of the Jeans length as 
M — > 0. Truelove et al. [106] have shown that unless the size of the grid, Ax, 
remains less than about 0.25 of the local Aj, the system is subject to spurious 
local fragmentation. 



3.2 Microscopic Parameters 

Complete expressions for plasma transport coefficients are given by |l3| . We are 
primarily concerned here with the magnetic and viscous diffusivities, which are 
conveniently given in terms of the electron and ion collision times r e and Ti 

_ 2.9 x IP' 2 T 3 / 2 1.7 A 1 ' 2 ^' 2 

Te ~ (A/10) n e Z S] Tl ~ (A/10) n e Z 3 (18) 

In eqn. (|lg|), T is given in degrees K, A is the Coulomb logarithm (A = 9.4 — 
1.151ogn e + 3.45 logT for T < 5.8 x lO 5 ^ and A = 15.9- 1.15 log n e +2.3 log T 
for T > 5.8 x 10 5 K), and A and Z are the ionic atomic number and charge, 
respectively. The electron density n e is expressed in cgs units; cm -3 . 

The magnetic diffusivity A^ introduced in eqn. (Q) can also be expressed in 
terms of the electron skin depth 5 e = c(m e /47rn e e 2 ) 1 / 2 = 5.4 x 10 5 n ( T 1 ^ 2 cm and 
T e as 

\ n = 6 f = 9.9 x 10 12 (yl/10) ^cm 2 s- 1 . (19) 
This numerical expression allows us to estimate the Lundquist number S 



5,^ = 0.0224^, (20) 
-V; nJ 2 (A/10) 

where the last expression holds in a hydrogen plasma. Under astrophysical con- 
ditions, S is enormous: for example, in ionized interstellar gas with T = 10 4 , 
n e = 1, B = 3 x 10" 6 G, and L B = 1 pc, S = 10 17 . With the exception of 
extremely dense, cold environments such as protostellar disks, or systems with 
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extremely small lengthscales, such as the outer layers of accreting neutron stars, 
Xq is always much less than the numerical diffusivity arising from discretization. 
The viscous diffusivity, or kinematic viscosity v, is 



kT 1 A x fO 8 T 5 / 2 
m, (A/10) n e A 1 / 2 Z 3 



r, = — — , , , _„ cm 2 s~ 1 . (21) 



The Reynolds number R with respect to an organized flow V on a lengthscale 
L is defined as LV/v, and is almost always extremely large, confirming our 
expectation that astrophysical flow can be structured over a wide range of scales. 

According to eqns. (|l^) and (^l|), the ratio of viscous to magnetic diffusivity, 
which is known as the magnetic Prandtl number P r , is 

v _ 1.4 x 1CT 5 T 4 
Pr= A^~ (A/10) 2 n e AW (22) 

Typically, P r 3> 1; in the example of fully ionized interstellar gas introduced 
above, P r — 10 11 . Therefore, the kinetic energy spectrum should be truncated 
on much larger scales than the magnetic fluctuation spectrum p5| , |66|p^ | The 
computations by Maron & Cowley |3(| show that the essentially constant shear 
of the velocity field on scales near the resistive scale generates tightly folded, or 
hairpin-like, magnetic structures. 

If the ion cyclotron frequency tu C i ~ 1.8 x lO^BZ/A and the ion collision 
time Tj satisfy the condition oj C iTi 3> 1, the viscosity is highly anisotropic with 
respect to the magnetic field. The viscous force acting on shear flow parallel to 
the magnetic field is reduced by a factor of (uj C iTi) 2 , while the parallel viscosity 
remains the same (see ]l3| for a full description). Using eqn. (f^), we see that 
UdTi is indeed large; in our numerical example, it is nearly 10 , implying very 
strong suppression of viscosity perpendicular to the magnetic field. Maron & 
Cowley |6(| have implemented the full tensor viscosity in their simulations of 
turbulent dynamos. 



3.3 Implications for Numerical Techniques 

In order for a numerical method to include diffusive processes, we would ideally 
require two things, namely (a) that it can resolve the physically important scales 
and (b) that it can handle the disparate time scales. Clearly, meeting each of 
these requirements alone is already a non-trivial task. 

Resolving the physically interesting scales means separating the intrinsic nu- 
merical diffusion existing in any scheme from its physical counterpart. Simulta- 
neously, it is desirable to cover scales much larger than the diffusive scales as 
well, in order to embed the problem in a realistic environment. Currently, studies 
are restricted to one of these regimes, for obvious reasons. 

There are various ways to parametrize the small-scale physics in the large 
scale approach. Shock capturing is a familiar example (for a detailed discussion 
see e.g. |5q| ). Many finite-difference codes apply a von- Neumann- Richtmyer ar- 
tificial viscosity |78|, spreading out the shock across several support points and 
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not affecting the solution in smooth regions. Artificial viscosity mimicks physi- 
cal viscosity in that it satisfies the Rankinc-Hugoniot conditions, but on a much 
larger scale, which renders the numerical scheme very robust at low computa- 
tional cost. It operates only where the flow is compressive. 

Galsgaard & Nordlund |36| implemented an artificial resistivity algorithm 
similar to artifical viscosity in their simulations of formation and dissipation of 
current sheets in the solar corona; see also Caunt & Korpi |^l[ . 

"Current sheet capturing" is intrinsically more complex than shock captur- 
ing. As we discuss in §4.3, the energy dissipated by a shock is independent of the 
magnitude of the diffusivity, but this does not hold for a current sheet. More- 
over, the topological evolution of the magnetic field is determined primarily by 
processes in current sheets, and this imposes additional requirements on how 
they are treated. Detailed, small scale studies could be useful for developing a 
parameterization of the effects of diffusion on energy balance and topology. This 
is reminiscent of the approach taken in large eddy simulations, which numeri- 
cally resolve the largest scales while treating the unresolved scales with a subgrid 
model (e.g. |l|§7|||). 

Godunov-type methods ]4l} | avoid broadening the shocks over many zones. 
They solve the Riemann problem defined by a set of advection equations for 
the physical variables and the physical states on both sides of the discontinuity. 
The upstream and downstream states are connected by a sequence of shocks or 
rarefaction waves, each of which satisfies the Rankine-Hugoniot conditions. As 
implementation of a full Riemann solver would involve determining the full wave 
structure and all wave propagation speeds, a variety of approximate Riemann 
solvers has been designed to simplify the process by disregarding certain wave 
types or linearizing the problem (see |38| for details). Riemann solvers utilize the 
wave nature of hyperbolic equations, however, the diffusive terms are parabolic, 
so that they cannot be included in the Riemann solver directly. 

An alternative approach has been taken in the so-called BGK-schemes J7], 
usin g the collisionless Boltzmann equation as a model of gas dynamics (see |8K ] , 
and Jl03 for MHD). 

Needless to say, the P r ^> 1 regime cannot be captured by numerical codes 
which solve the ideal fluid equations (77 = v — 0), since the numerical diffusivities 
for magnetic field and velocity are about the same size and occur at about the 
same scale, i.e. the grid scale. This is an argument for including viscosity and 
resistivity in the equations, in which case the physical diffusivities must be larger 
than their numerical counterparts, and their scales must be well separated. This 
is the approach taken by Maron & Cowley J6(| , who used a spectral code to study 
turbulence with widely separated, but resolved, viscous and resistive scales. 

Several efforts have been made to connect the spatial scales by adaptive 
mesh refinement techniques ||,|J for large eddy simulations p7|] , general hyper- 
bolic systems [|J , incompressible non-ideal MHD (34) and for compressible MHD 
[110 i|. However, turbulent flows by definition connect scales over the whole do- 
main, which complicates finding a suitable refinement criterion and may in the 
end yield only modest savings of time. 
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We have seen in §3.1 that timescales in astrophysical problems can differ by 
several orders of magnitude. Including diffusive processes aggravates the prob- 
lem, as their time scales are generally orders of magnitude longer than the dy- 
namical times. Any disparity in timescales leads to a stiff problem. Explicit 
methods - often the first choice because of their low time and memory needs - 
advance the solution by a time step At which is given by the so-called Courant 
- Friedrichs - Levy (CFL) condition p8| 

At < — (23) 
c 

with the distance between two support points Ax and a characteristic prop- 
agation speed c. Disparate characteristic speeds then lead to severe time step 
restrictions. Moreover, the slowly varying components may introduce numerical 
errors fsofl . 

One solution to the problem is to choose a short timestep for the fast pro- 
cesses, and update the variables controlled by slow processes less frequently. 
Mac Low et al. used this so-called subcycling in their implementation of 
ambipolar diffusion in the ZEUS code. Depending on the problem, the critical 
processes often can be broken out and treated implicitly, while the remaining 
equations can be treated explicitly Q. Fully implicit schemes (e.g. can be 
unconditionally stable and independent of the time step chosen. Their main limi- 
tation lies in the tremendous computational needs in the multi-dimensional case. 
For a discussion of implicit methods and for a comparison between explicit and 
implicit shock capturing methods see 50 10£]. For a discussion of higher order 



finite difference schemes and their application to MHD turbulence simulations 
see 0. 

As this brief discussion shows, a great variety of techniques can be brought 
to bear on MHD turbulence problems. Which technique is best depends very 
much on the problem at hand. 



4 Results from Theories and Simulations 

Theory offers an interpretive framework for simulations, while numerical ex- 
periments test turbulence theory, and can inspire new developments. Here, we 
discuss amplification of a field by turbulence in the context of dynamo theory, 
magnetic fluctuations, the formation of current sheets, and the dynamical effects 
of turbulence on large scales. The first three of these problems are closely related. 
The fourth is topical in view of recent simulations of star forming regions. 

4.1 Turbulent Amplification of a Weak Field 

Amplification of a weak magnetic field by turbulence is one of the main compo- 
nents of dynamo theory. A successful astrophysical dynamo theory must demon- 
strate that a dominant, large scale magnetic field can be generated by small 
scale turbulence, on a timescale which is virtually independent of the Lundquist 
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number S in the limit S — > oo. This is a highly nonlinear problem, which involves 
multiple scales. It probably cannot be solved without numerical simulations, but 
in view of the extreme parameters involved it is unlikely to yield to brute force 
alone. 

The growth rate of magnetic energy in a periodic domain is given by eqn. 
( p"2| ) . Although the resistive term appears as a sink, some resistivity is necessary 
to effect irreversible topological change. Amplification occurs through work done 
by the flow. Ingeneral, there is also a surface integral — c J dS ■ E x B on the 
RHS of eqn. (|l2T), which represents an energy flux through the boundaries. 

Equation ( |12|) says nothing about the structure of the field, and, in particular, 
whether magnetic energy is concentrated at large or small scales. If we formally 
set \q = 0, the ratio of fieldstrength to line length is preserved by incompressible 
motions. The only way to lengthen fieldlines in a finite volume is to wind or 
tangle them. Therefore, we expect that most of the energy will be in the rms 
field, rather than in the mean field. For example, if the protogalactic magnetic 
field were 10~ 17 G [flO), the fieldlines would have been lengthened by more than 
a factor of 10 11 in the course of amplifying the field to its present strength. 
Yet, the large scale and small scale components of the galactic magnetic field 



are observed to be roughly equal [ 1 13 1 . It is not enough to argue that magnetic 
forces will prevent the field from becoming tangled on small scales. If the lines 
aren't somehow lengthened, the field won't be strengthened. 

Of course, we are interested in — > 0, not = 0. The fieldlines are allowed 
to break, and the question is whether they do so in a way which prevents power 
from piling up at the resistive scale and peaking instead at the large scale. 
Because S is so large, topological constraints are strong, and we expect resistive 
boundary layers to form. 

Let's see how theories and simulations of dynamo theory address this prob- 
lem. The most influential and detailed dynamo theory is the so-called mean 
field theory, which was proposed by Parker p3| , and extended and formalized 
by Steenbeck, Krause & Radler |97|) . Mean field theory is based on the idea 
that a large scale magnetic field (B) can be generated from the average induc- 
tive electric field (E) associated with small scale, helical velocity and magnetic 
fluctuations Sv, SB. The dynamo property of small scale, helical fluctuations is 
known as the a effect. There is also diffusive transport of (B), which is known as 
the (3 effect. Although a and (3 are in general tensors, it suffices for our purposes 
to write them as scalars. The induction equation for the mean field is 

& = -cVx{E), (24) 

where 

(cE) + (V)x(B) = \ n Vx(B)-(5vx8B) = a(B) + (\ n + (3)Vx(B)+.... (25) 



In eqn. (|25|), (V) is the large scale velocity field, and the represent additional 
terms involving moments of the turbulent fluctuations and successively higher 
derivatives of (B) (see Moffatt Note the correspondence between eqns. 

(HJ) and (||) and eqns. (0) and (|)^ 



Simulations of Magnetic Fields 



13 



If a and /3 are independent of (£?), eqn ( p4| ) is linear in (B). Physically, 
this corresponds to neglecting the back reaction of the magnetic field on the 
turbulent velocity, so the linear case is sometimes called the kinematic case. 

The linear version of eqn. ( p4| ) has solutions which vary exponentially in 
time. At large shear and low wavenumber k, the growth rate is approximately 
the geometric mean of the shear rate V and the turbulent frequency ka. 

Standard mean field dynamo theory does not directly address the evolution 
of the small scale field, and the role of resistivity, on which the fate of the small 
scale field depends, is left somewhat implicit. Usually is dropped from eqn. 
(p4|), since f3 is assumed to be much larger. A calculation by Moffatt of the 
a effect due to helical Alfven waves demonstrates that a = unless there is a 
phase difference between Sv and SB. In Moffatt's model, the phase lag arises 
from magnetic diffusivity. At small A^, Q oc Aj7, implying slow dynamo action 
at large R m . 

In more general calculations of a, the diffusive damping time (fc 2 A^)^ 1 is 
replaced by the correlation time r of the turbulence. At large R m , this is much 
shorter than the resistive time, and implies large growth rates for the dynamo. 
But (3 does not distinguish between converting large scale field to small scale 
field and actually destroying field, so the behavior of the small scale fields are 
obscured by this treatment. Kulsrud & Anderson [^5| used the tools of mean 
field theory itself to follow the growth of the small scale field. They showed that 
at large R m , the rms field (B 2 ) 1 / 2 grows much faster than the mean field (B), 
and is dominated by small scale fluctuations with a fc 3 / 2 power spectrum. 

A completely different perspective on kinematic dynamo theory emerges from 
solutions of the full magnetic induction equation for prescribed, chaotic flow. 
These studies, which are fully described up to 1995 in the volume by Childress 
& Gilbert p3| , use a combination of numerical simulations at large R m and 
analytical theory to establish connections between the properties of the flow 
and the properties of the magnetic field. Although there are known examples of 
flows which amplify the magnetic field at a rate independent of R rn , as originally 
shown by Galloway & Proctor [B5J, the fields produced by these so-called fast 
dynamos are highly intermittent in space and fluctuate in sign, confirming at 
least qualitatively that small scale fields are dominant in kinematic dynamos at 
high Rm. 

The kinematic theories are modified by dynamical effects. It was expected 
on general grounds that the dynamo would saturate as the magnetic field ap- 
proached equipartition with the kinetic energy. A saturation process known as 
"Alfvenization" was first described by (Pouquet et al |8g|| ). Saturation occurs 
because the small scale fluctuations increasingly resemble Alfven waves as the 
mean field grows. In an ideal medium SB — > ±Sv in Alfven units, so their cross 
product also tends to zero, eliminating the a effect (also shown in the calcula- 
tion by Moffatt). Alfvenization was first identified by Pouquet et al. (88) using 
a spectral closure scheme, and later derived using quasilinear theory by Gruzi- 
nov & Diamond [Q. In the model of Pouquet et al., kinetic helicity injected 
at the forcing scale generates magnetic helicity, which undergoes both a direct 
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cascade, to high k, where it is resistively dissipated, and an inverse cascade, to 
low k. Alfvcnization saturates the dynamo once the large scale field is roughly 
in equipartition with the kinetic energy. In the model of Gruzinov & Diamond, 
saturation by Alfvenization takes place when the ratio of kinetic to large scale 
magnetic energy is 0(R m ), confirming the dominance of small scale fields at 
large R m . The difference between these calculations may lie in their treatments 
of diffusion. Gruzinov & Diamond explicitly include \q, while Pouquet et al. 
took all the diffusivities to be eddy diffusivities, which are much larger than 
molecular diffusivities. 

The role of kinetic helicity injection has been tested in simulations. Compu- 
tations by Maron & Cowley [ |66| , without kinetic helicity injection, show a peak 
in spectral power at the resistive scale, as predicted by Kulsrud & Anderson ]55|) . 
Simulations by Maron & Blackman Q show that a large scale field is generated 
when kinetic helicity is injected, but that the growth rate is proportional to the 
resistive time. Nonlinear simulations with helical forcing by Brummell et al. |]l7f 

— 1/2 

show that the magnetic lcngthscale is proportional to R m , while kinematic 
investigations testing the sensitivity of fast dynamo action to kinetic helicity 
have not found an inverse cascade Q . 

Thus, we see that much hinges on the resistive time, or, more broadly, the 
mechanism by which small scale fields are dissipated. Turbulent diffusion of (B), 
parameterized by the (3 effect, merely stands for spectral transfer out of the large 
scale. The possibility remains that turbulence could mix the field to the small 
scale, efficiently destroying it. However, numerical models by Cattaneo et al. 
[ pp| , and subsequent analytical work by Kim |32| suggest that the mixing is self 
limiting due to Lorentz forces. 

The possibility that the small scale fields might escape through an open 
boundary rather than being destroyed in situ, leading to growth of the large 
scale field on the convective timescale rather than the resistive timescale, was 
raised by Parker |56| and Blackman & Field |ll| Up to now, neither analytical 
studies |)2j nor numerical models jl5| have confirmed this idea. 

Astrophysical dynamo theory is at a critical juncture. The fundamental 
premise of mean field dynamo theory - that kinetic helicity injected at small 
scales can drive the growth of a magnetic field at large scales - is supported 
by a variety of analytical treatments and numerical investigations. Mean field 
theory is a felicitous outcome of turbulence theory in the sense that it can be 
applied just by calculating a few parameters (or functions): a, /3, and the large 
scale shear. However, neither theory nor simulations have fully come to grips 
with the large value of R m and its implications for the small scale fields. From 
a numerical point of view, computations at large R m which can capture current 
sheets as well as accommodating a variety of boundary conditions would seem 
to be a prerequisite for either validating or superseding mean field theory. 

4.2 Magnetic Fluctuations 

Steady hydrodynamic turbulence, uninfluenced by boundaries, has been charac- 
terized in two regimes. Subsonic turbulence is incompressible, and follows the 
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Kolmogorov scaling (E(k) oc fc~ 5 / 3 ). Supersonic hydrodynamic turbulence is 
compressible, and better described as Burgers turbulence (E(k) oc fc~ 2 ). (These 
power law spectra, and all subsequent ones, describe only the so-called inertial 
range, in which there is no driving or microscopic dissipation, only nonlinear 
spectral transfer.) 

Both these limiting forms of turbulence are spatially intermittent in the sense 
that the variances of quantities such as the rate of energy dissipation undergo 
substantial fluctuations from point to point. In Burgers turbulence the basic 
structures are shock waves, and the large fluctuations and gradients occur prin- 
cipally in the shock fronts. In Kolmogorov turbulence, the basic units are ed- 
dies, and intermittency is represented by the concentration of vorticity into small 
structures. Intermittency does not affect the energy spectrum of Burgers turbu- 
lence (the Fourier spectrum of a shock itself is fc -2 ), but is thought to steepen 
the Kolmogorov spectrum. This occurs through the local enhancement of the 
dissipation rate at small scales fl57fl . 

Hydrodynamic turbulence in unstratified, nonrotating systems is isotropic. 
A uniform component of magnetic field causes turbulence to be anisotropic. 

An ideal, adiabatic fluid with a uniform component of magnetic field has 
three linear modes: the shear Alfven mode, and the fast and slow magnetosonic 
modes. The velocity field Sv of the Alfven mode, is solenoidal, and this mode is 
driven purely by magnetic tension. The dispersion relation is u> = k\\VA = u>a, 
where fc|| is the component of k parallel to the large scale magnetic field. The 
fast and slow magnetosonic modes are compressive, and driven by a combination 
of gas pressure, magnetic pressure, and magnetic tension. 

Since the shear Alfven mode is noncompressive, it is reasonable to ask whether 
a strong (/3 < 1) magnetic field permits a regime of pure shear, supersonic tur- 
bulence. The answer is no. 

First, suppose k± = 0. Then, the nonlinear magnetic pressure gradient associ- 
ated with the fluctuating transverse magnetic field drives a compressive parallel 
flow, causing an Alfven wavetrain to steepen nonlinearly into a train of weak 
shocks [p6| , much as a sound wave steepens. The only exception is the infinitely 
long, circularly polarized wave, which is an exact solution of the ideal MHD 
equations. The outcome of the evolution of an ensemble of parallel propagating 
Alfven waves is a series of shocks, similar to Burgers turbulence, with k~ 2 kinetic 
and magnetic energy spectra |]39| . The steepening time depends on amplitude 
as (Sv/va) _2 , while the steepening time of an acoustic wave is proportional to 
(Sv/vs) -1 ■ In this sense, shear waves survive longer in a low (3 plasma than 
compressive waves of the same amplitude, by a factor of order P~ 1 / 2 (va/5v), 
but the outcome in either case is a compressive flow dominated by shocks, and 
the difference in timescales is only significant for small wave amplitude. 

Although Burgers turbulence was originally based on ID flow, 3D simulations 
of supersonic MHD turbulence show similar rapid evolution to a shock dominated 
flow. The distribution of shock strengths in simulations of driven and decaying 
turbulence, with and without magnetic fields, has been studied by Smith et al. 

mi- 



16 



Ellen G. Zweibel et al. 



The steepening of Alfven waves, which can be viewed as a cascade in fcy, is 
suppressed if (3 > 1 . Under the assumption of incompressibility, and in contrast 
to the highly compressible low (3 case, interactions occur only between oppositely 
directed wave packets. In particular there is no self interaction of the kind which 
leads to steepening. 

The incompressible shear Alfven wave cascade is highly anisotropic, with 
power transferred much more rapidly in k± than in fen . In physical space, the 
correlation length transverse to the mean field is much shorter than the cor- 
relation length alon g it. These anisotropic states can be described by so-called 
reduced MHD [ 101 1 , which is a quasi-2D approximation describing elongated 
magnetic structures. It is often applied to (and, indeed, was derived for), low 
(3 plasmas, but it neglects the parallel steepening which leads to shocks. This, 
as we said above, is a good approximation for long parallel Alfven transit times 
and small turbulent amplitudes. 

There is not yet complete consensus on the cascade itself. Weak MHD tur- 
bulence, in which an individual Alfven wave packet survives for many periods, 
can be viewed as weakly perturbative resonant interactions between multiple 
waves. Sridhar & Goldreich |9(| argued that the dominant interaction is a four 
wave interaction, which results in a spectrum E{k) oc fc -7 / 3 . Others |7j],|||38| 
claim that the dominant interactions are three wave interactions, and predict 
E(k) oc k~ 2 . 

In strong MHD turbulence, wave packets survive for only about one wave 
period. Goldreich & Sridhar developed the concept of a critically balanced 
cascade in which the wave frequency to a is the same as the nonlinear frequency 
k±v±. This, together with the requirement of constant spectral energy flux ar- 
gument k±v^_, leads to a Kolmogorov spectrum E{k) oc fc~ 5 / 3 and wavenumber 

— 1/3 

anisotropy k\\/k±^ oc k ± . Other theories of strong MHD turbulence predict 
E(k) oc fc~ 3 / 2 , including the original isotropic Iroshnikov-Kraichnan theory fl% . 

Numerical simulations, if free of confounding computational effects, could 
play a role in resolving these disagreements. Biskamp & Miiller [[[o) and Cho 
& Vishniac |25|| find E{k) oc fc~ 5 / 3 over about one decade in k space. Maron 
& Goldreich |67| slightly extended the inertial range and found E{k) oc fc~ 3 / 2 . 
These authors argue that the spectrum is flattened because of intermittency. 
The crucial difference between the HD and MHD cases is that energy cascades in 
MHD only when oppositely directed wave packets collide; the small filling factor 
associated with intermittency reduces the collision rate, more than offsetting the 
enhanced dissipation associated with small structures. 

Larger computations, with an extended inertial range, should shortly become 
available. The differences in the spectra obtained with different codes and under 
different forms of driving (Biskamp & Miiller studied decaying or isotropically 
forced turbulence, Cho & Vishniac studied isotropically forced turbulence, and 
Maron & Goldreich studied anisotropically forced turbulence) are at this point 
comparable to the theoretical disagreements. A joint exercise in which differ- 
ent groups simulated identical problems and implemented identical diagnostics 
might be quite enlightening. 
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A variety of processes can terminate the cascade at short wavelengths (see |39[ 
for a recent discussion). These include fluid effects, such as ion-neutral friction, 
and kinetic effects such as ion gyroresonance absorption and electron Landau 
damping j9l) . The extension of the magnetic spectrum to scales below the viscous 
cutoff in a high Prandtl number plasma is discussed by |66) and [^4| . 



4.3 Intermittency: Current Sheets and Flux Tubes 

Simulations of MHD turbulence show strong intermittency in the distribution 
of magnetic field gradients, or current. In the weak field, high Prandtl number 
simulations of |6q] , the fluctuating field is tightly folded. These results suggest 
that we look at current sheets. 

Current sheets and filaments are sites of intense local heating, and possibly 
particle acceleration, which makes them interesting in their own right. It is quite 
likely that in most astrophysical systems, significant departures from flux freez- 
ing occur only in these regions, so their existence is important in the evolution 



of magnetic field topology, and for dynamo action, (see §4.1). Reconnection at 
magnetic X-points drives strong, small scale jets. Simulations by Galsgaard & 
Nordlund |37j show highly time variable dissipation rates associated with the 
formation and disruption of large scale current sheets. 

Current sheets are a prerequisite of, and accompany, magnetic reconnection. 
Reconnection is a resistive process which occurs on a timescale intermediate be- 
tween the Ohmic time and the AlfVen time. The reconnection time r rec is often 
parameterized by the Lundquist number S; r rec scales as S p with < p < 1 . If 
p = 0, the reconnection is said to be fast (in correspondence with the definition 
of a fast dynamo). In Sweet-Parker reconnection, p = 1/2 |p5|| , and resistive 
tearing modes in slabs and cylinders have p — 3/5, and p — 1/3, respectively ||. 
The formation and behavior of current sheets has been studied for many years. 
Most of this work has focussed on plasmas in or near magnetostatic equilibrium, 
which can be expected if /3 <C 1 and the plasma is either allowed to relax or 
is driven at frequencies far below its characteristic AlfVen frequency. Such con- 
ditions hold, for example, in stellar coronae, in which the magnetic fieldlines 
are tied to the underlying photosphere, and evolve quasistatically in response to 
slow photospheric motion. 

Current sheets form as a plasma seeks equilibrium while obeying the topolog- 
ical constraints imposed by flux freezing |75|] . Topological constraints can arise 



from features such as null points [ 102 1 or separatrices [[32 115| ,60|- Sweet [ 103 1 



and Parker |8j] suggested that shearing and winding the footpoints of an initially 
uniform field will create current sheets. Even a simple sinusoidal deformation of 
the boundaries of a sheared magnetic slab can induce the formation of a current 
sheet (||. 

As a simple example of how current sheets can form, consider the magnetic 
field of a pair of neighboring, uniformly magnetized superconducting spheres, 
with antiparallel magnetic dipole moments, surrounded by a perfectly conduct- 
ing, zero pressure plasma. Some fieldlines have both endpoints on the same 
sphere, and others have one endpoint on each sphere. 
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Now suppose that one sphere rotates on its magnetic axis by an angle <f>. 
The fieldlines which are not connected to the other sphere also rotate by cf>, but 
the fieldlines which are connected to the other sphere cannot rotate uniformly 
because one endpoint remains fixed. The separatrix surface which divides the 
two domains of magnetic connectivity becomes a current sheet: the field outside 
it is sheared while the field inside it is not. 

There is still no comprehensive theory for exactly when current sheets form, 
and how their properti es re flect the underlying magnetic topology. Nevertheless, 



analytical arguments [ 108 ,61| and numerical experiments |73| , |6l| , |36| show that 
shear, or equivalently current density, grows exponentially rapidly in line tied 
magnetic fields with randomly moving footpoints. This has some of the same 
consequences as an exact singularity, which cannot form in any case in even a 
slightly resistive plasma. 

These results are not fully applicable to MHD turbulence, in which the mag- 
netic field is not in equilibrium. Nevertheless, current sheets form in flows with 
magnetic null points [jS7| . The underlying mechanism for the growth of the cur- 
rent is the random stretching of neighboring fieldlines, which leads to large cross- 
field shear. 

What are the consequences of current sheet formation for dissipation of mag- 
netic energy? That is, if energy per unit volume e is being added to a system by 
a driving process, under what circumstances can this energy can be dissipated 
in current sheets? 

Consider a current sheet of area L 2 and width S. Assume the current sheet 
structure is given by the Sweet-Parker theory, so that 6 = LS^ 1 / 2 = {LAq/va) 1 ^ 2 ■ 
If N(L)dL is the number density of current sheets with L between L and L + dL, 
then the rate at which energy is dissipated by current sheets in this size range is 

dw = —^±N(L)L 2 5dL = —N(L) (L 3 \ n v A ) 1/2 dL. (26) 

47T L All 

The rate of energy dissipation by all current sheets is obtained by integrating 
eqn. ( p6| ) over L: 

* = J ( A ^) V2 / dLN(L)L^ 2 . (27) 

If e is independent of Aq as A^ — > 0, then J dLN(L)L 3 / 2 must scale as A n 1/>2 in 
this limit. This would imply an unbounded number of current sheets as An 0, 
which is virtually equivalent to a pileup of magnetic energy in fluctuations at 
the resistive scale. [] 

There are alternative possibilities. One is that the current sheets are not 
structured according to the Sweet-Parker picture. If S scaled as A^, as it does 
(up to a logarithmic factor) in Petschek's model of fast reconnection (reviewed in 

1 If we make a similar argument for the efficiency of energy dissipation in shocks, the 
result is independent of viscosity, because the shock thickness is directly proportional 
to viscosity. Therefore, the inviscid limit does not require an infinite number of shocks 
to dissipate the energy input by driving. 
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p5|), then J dLN(L)L 3 / 2 could be independent of A^- Note, however, that while 
in the Sweet-Parker theory, magnetic energy is converted to kinetic energy and 
to heat at about the same rate, in Petschek's theory magnetic energy is converted 
predominantly to kinetic energy, and that the realizability of Petschek's model 
has recently been questioned |9,107|. 

The second possibility is that the magnetic and velocity fields adjust them- 
selves such that the driving just balances the damping, as they do in mechanical 
systems such as linear harmonic oscillators. Thus, the energy input rate e would 
depend on Xq, while the applied force itself, presumably would not. 

Dedicated simulations of current sheets and reconnection regions themselves 
allow us to study the small scale aspects of the problem. It is equally important 
to understand how current sheets are embedded in the overall flow. This has 
already been done to some extent for laboratory experiments [107]. Ultimately, 
it may be possible to parameterize magnetic reconnection in simulations on larger 
scales. 

Finally, we mention an even more extreme form of intcrmittency: the seg- 
regation of the magnetic field into thin tubes. This is observed in the solar 
photosphere, where magnetic flux tubes, or sheets, just a few hundred km thick 
are seen at the borders of convective cells. Flux tubes can form only in high (3 
plasma, so that gas pressure can balance magnetic pressure at the tube walls. 

Early attempts to explain photospheric flux tubes relied on a two stage pro- 
cess, proposed by Parker, for concentrating a diffuse field into organized struc- 
tures. According to this scenario, thermal convection would sweep the field to the 
borders of convective cells, concentrating the field up to equipartition with the 
flow. Then reduced convective heat transport would cool the gas, "collapsing" 
the tube and bringing it to equipartition with the thermal plasma. 

An alternative possibility is that intermittent fields are generated in situ. 
It has been shown that small scale turbulent thermal convection operates as a 
dynamo, and generates a magnetic field with a broad tail of high energy features 

m 

It is unclear whether a dynamo operating in a high (3 plasma should always 
produce flux tubes. In the solar case, the turbulence is organized by stratification 
and the thermal gradient, and is subsonic. If the turbulence were supersonic and 
the field were amplified to equipartition, then flux tubes could not persist. 



4.4 Dynamical Effects of Turbulence 

The stresses associated with turbulence can exert dynamical forces. In MHD tur- 
bulence, these forces include both magnetic stresses and the Reynolds stress. Ob- 
servations suggest that these stresses are important in the interstellar medium, 
and especially in molecular clouds. 

It was realized early that supersonic motions are nearly ubiquitous in molec- 
ular gas, and pose a problem. If the motions are associated with gravitational 



collapse, the implied star formation rate would be very high [ 1 1 1 L but if the 
motions are turbulent, they should be dissipated very quickly §42|| . This co- 
nundrum, together with the expectation that reasonably strong magnetic fields 
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should be present, prompted analytical studies of Alfvenic turbulence under 
molecular cloud conditions |p], |l 14 32 . This work demonstrated that even when 
Ma <C M and the waves are purely transverse, ion-neutral friction and nonlinear 
steepening are strong damping mechanisms, which limit the lifetimes of Alfven 
waves to a few 10 6 yr or less under molecular cloud conditions. 

Even so, under the assumption that the waves might be replenished by 
some energy source, various aspects of the theory of self gravitating clouds sup- 
ported by Alfven wave stresses were developed and applied to model clouds 
[ po[|3^ , |69| , |68|j7C|l . According to this theory, which is based on the the weak turbu- 
lence approach pioneered by Dewar [p9[ , Alfven waves have an isotropic stress 
tensor P = P W I. The waves pressure - density relation is P w oc p 1 / 2 in a strat- 
ified, static medium, while P w oc p 3 ^ 2 under slow, spatially uniform changes in 
density with time. The negative polytropic index (n = — 1/2) of the P— p relation 
leads to a large center to surface pressure contrast, in accord with observations 
(sec [[70| for a good discussion). 

The theory of wave supported clouds is most relevant to objects which are 
slightly magnetically supercritical, that is, the ratio of mass to magnetic flux is 
slightly too large for magnetostatic support. If the clouds were subcritical they 
would be supported by the DC magnetic field. If they were highly supercritical, 
the amplitude of the turbulence required to sustain them would be so large that 
the weak turbulence theory would probably not apply, although they could still 
be turbulcntly supported. ^From an observational viewpoint, the most likely 
candidates for wave support are dense cores with sonic or mildly supersonic 
linewidths Q. 

The picture of strongly supersonic turbulence which emerges from numerical 
simulations of molecular clouds is quite different from this analytical picture, pri- 
marily because, as we discussed in §4.2, generation of compressive disturbances 
is entirely unavoidable. These compressive disturbances, or shocks, sweep up 
most of the mass into thin layers in which the local Jeans length is small. These 
layers, and the even denser structures formed where layers intersect, collapse if 
the domain is magnetically supercritical, whether the turbulence is freely decay- 
ing or maintained in a steady state by driving J471 , ^l| . The driven models are 
globally stable in the sense that Aj computed using the total turbulent energy 
and mean density is larger than the length of the domain. The small, high den- 
sity structures collapse because the turbulent power at wavelengths less than the 
local Jeans length is relatively small (see eqn. (|l7|)), and because the intensity of 
turbulence inside the cores is comparable with the density outside. Nevertheless, 
the bound structures are small, and the role of numerical diffusion within them, 
which damps the turbulence and removes the support by the DC magnetic field, 
has not been quantified. Therefore, while the formation of turbulently supported 
clumps has not been found in the numerical models, their existence cannot be 
ruled out based on the present models. 
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5 Ambipolar Drift 

As we have discussed in the preceding sections, both theoretical and computa- 
tional aspects of astrophysical magnetohydrodynamics are dogged persistently 
by the enormous value of R m . There is one situation, however, in which the mag- 
netic diffusivity is large whether or not it is enhanced by turbulence: in weakly 
ionized gas, the magnetic field and plasma drift with respect to the neutrals. 
This so-called ambipolar drift is not entirely equivalent to resistive diffusion, 
because it preserves magnetic topology (the field is frozen to the plasma) , but it 
does change the ratio of magnetic flux to mass, as originally pointed out in |72| , 
and dissipates fluctuations on small scales. 

Full treatment of ambipolar drift requires solving the equations of MHD 
for the plasma and neutrals treated as distinct species coupled by collisions, 
ionization, and recombination pc| |. At the low ionization fractions expected in 
molecular clouds, the ion Alfven speed va% = B / ' ^JAitpi can be several orders of 
magnitude larger than the other characteristic speeds in the problem, reducing 
the maximum possible timestep by a similar factor. This and other numerical is- 
sues related to implementation of ambipolar drift are discussed by |33, 105 64| { 



On timescales longer than the ion - neutral collision time r,- n , and at low 
ionized mass fractions, the ion-neutral drift Vd is determined by balancing the 
Lorentz force on the ions against the frictional drag by neutrals 

(V x B) x B _ v 2 M r m f 
v d = = — /, \2x>) 

VKPiVin L B 

where / is a unit vector in the direction of the Lorentz force. The second equality 
in eqn. ( p8j ) can be taken as a definition of the magnetic lengthscale Lb- The 
quantity uJyTj n has units of diffusivity, and from now on we refer to it as Xad- 
With this definition, vd = Xad /Lb- Note that pi/ri n = p n /T n i and p n » p, so 
v M T in — v A T n i, and Xad is often expressed in the latter terms. 

The ambipolar Reynolds number Rad is defined as the ratio of the bulk 
velocity v to the drift velocity Vd 

R AD = ^i, (29) 

*AD 

and measures how well the magnetic field is frozen to the bulk fluid on scale L B 
[ |ll2| . Setting Rad = 1 correctly predicts the thickness of shocks in which the 
main dissipation mechanism is ion- neutral friction |3l|] , the wavelength at which 
Alfven waves are critically damped jpq] , and the minimum size of an eddy which 



can wind up a magnetic field [112] 



Simulations of supersonic MHD turbulence [|82| with varying strengths of 
ambipolar drift shows that, on average, the drift smooths out the current, and 
reduces the rms Lorentz force. 

We expect (with one caveat discussed at the end of this section) that the 
magnetic field should show very little structure below the scale L min = \ad/v 
at which Rad = L Consider a volume L 3 of turbulent gas with mean density 
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p, turbulent velocity v, and mean magnetic field B. If we introduce the crossing 
time Td = L/v and the Alfven Mach number Ma = v/va, then we can write 

~ir = w A — d - m 

Expressing v in units of km s _1 as V5, and L in pc as L pc , t<i~3x 10 13 L pc /v5 
s. For the particle species expected in molecular clouds, t„, ~ 5 x 10 8 ni X s. With 
the ionization equilibrium relation rn ~ lO -5 ™; 1 / 2 (cgs; (7l)), eqn. (|3fj| ) becomes 

Comparing eqn. ( |3l| ) with the grid spacing Ax/L = 1/N, we see that a moder- 
ately large numerical simulation (say 256 3 ), which is not highly super- Alfvenic 
or extremely dense, should resolve almost all the magnetic structure associated 
with turbulence. 

Equation ( |30| ) can also be written in terms of the critical magnetic field 
B c = 2nG 1 / 2 pL and the free fall time r ff = {AnGp)- 1 / 2 - 5.5 x 10 14 n^ 1/2 s as 



Bj 4r 2 f \Bj Ar ff ' 



(32) 



where the last step holds for virial equilibrium; t// ~ Td- Of course, eqn. ( po| ) is 
entirely independent of G. 

The modest value of Rad in molecular clouds has been commented upon 
and exploited by Hg). 

Ambipolar drift is essentially diffusive when the magnetic field has a well 
ordered, well combed structure. Because the ambipolar diffusivity is proportional 
to B 2 , the diffusion is nonlinear. Sharp fronts can be generated in the vicinity 
of magnetic nulls, and minima are steepened pp|,|ll|,pl. 



6 Summary and Future Agenda 

The theme of this paper is how numerical simulations can best cope with the 
enormous magnetic Reynolds numbers encountered in most astrophysical prob- 
lems. We argued that, because R m (and the Lundquist number S) are so large, 
resistive effects occur primarily in thin sheets or filaments, which must be treated 
accurately in order to follow the topological evolution of the magnetic field. 

In §2, we reviewed the magnetic induction equation, giving the solution in 
the ideal limit, and derived equations for magnetic helicity and magnetic energy. 
We tested helicity conservation in the ZEUS code by simulating the evolution of 
the ideal kink mode, and studied fluctuations of helicity in a simulation of a tur- 
bulent cloud. The first of these problems is especially suitable for benchmarking 
numerical codes. That is, it yields a quantitative estimate of the numerical resis- 
tivity in the code, from which one can compute the resistive timescale (at least 
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in smooth regions, on large scales) and assess whether it is much longer than 
other timescales of interest. One can also compare the performance of different 
codes. 

In §3, we discussed some important dimensionless parameters in astrophys- 
ical MHD problems. Wc verified the large size of the Lundquist number S and 
magnetic Reynolds number R rn under most conditions, and showed that the mag- 
netic Prandtl number, the ratio of viscous to magnetic diffusivity, is also usually 
large. This implies that magnetic structure will generally extend to smaller scales 
than velocity structure. We briefly summarized a variety of numerical techniques 
from the perspective of the large S limit. 

In §4, we discussed theoretical and numerical results on several basic prob- 
lems. In §4.1 we discussed the amplification of magnetic fields by turbulence, and 
argued that the successful operation of a large scale dynamo depends on fast dif- 
fusion of the magnetic field. In §4.2 wc discussed MHD turbulence itself. In §4.3 
we discussed current sheets, which may provide the fast diffusion necessary for 
a dynamo, and may be produced in turbulence as a form of intermittency. In 
§4.4, we discussed the dynamical effects of turbulence, and its importance in 
molecular clouds. 

In §5, we discussed ambipolar drift as a mechanism for increasing the mag- 
netic diffusivity in low density, weakly ionized gas. It is now possible to simulate 
turbulent molecular clouds at a resolution which accounts for all magnetic struc- 
ture down to the ambipolar diffusion scale. Under certain conditions, however, 
ambipolar drift can mediate the formation of current sheets, possibly leading to 
rapid magnetic reconnection. 

It is fortunate that a number of different groups are simulating MHD tur- 
bulence under astrophysical conditions, using a variety of codes and techniques. 
There is substantial agreement on a number of issues, such as the short decay 
time of supersonic turbulence, but disagreement on others, such as the spectrum 
of strong MHD turbulence. A dedicated effort at benchmarking, in which differ- 
ent groups simulated identical problems and implemented identical diagnostics, 
would provide some perspective. We suggested in §4.2 that such an exercise be 
carried out for strong MHD turbulence. 

It is unlikely that it will ever be possible to adequately resolve current sheets 
and global dynamics in a single simulation. Local studies of current sheets and 
magnetic reconnection are essential in capturing the physics of these layers. It is 
equally important to understand how they are embedded in the overall flow, as 
this constrains their properties, sets the boundary conditions for reconnection, 
and is necessary for understanding how to parameterize the effects of current 
sheets in global simulations. The inability to deal adequately with this small 
scale structure is the greatest present challenge to understanding magnetic field 
evolution in turbulent flows. 
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